#include "VectorFunction.H"
#include "ExampleFunctions.H"

using namespace std;

int main(int argc, char *argv[])
{
    TestExactSol2D TS2D;
    const int NumComponents = TS2D.getNumComponents();
    IntVect lo(IntVect::Zero);
    RealVect origin(IntVect::Zero); 
    const int numLevel = 5;
    clock_t start, end;
    cout << "Served for 2D situation. "<<endl;
    cout<<"1. CPU Time Test for Setting Values on Boxes"<<endl;
    for (int i = 0 ; i < numLevel ; i++)
    {   
        ///Pick different origins for different boxes.
        RealVect origin_new(i,-i);
        IntVect hi(IntVect::Unit*(pow(2,i+7) - 1));
        Box B(lo,hi);
        FArrayBox FB(B,NumComponents);
        Real dx(pow(2,-i-4));
        start = clock();
        TS2D.setFAB(FB,origin_new,dx);
        end = clock();
        cout << "For sizes " <<  pow(2,i+7) << "*" << pow(2,i+7) << ", CPU time is "
             << (double)(end - start) / CLOCKS_PER_SEC<< " s." <<endl;
    }
    cout << " " << endl;
    cout << "2. CPU Time Test for 4th Discrete Laplacian: " << endl;
    Vector<Real> errU, errV, errP;    
    Vector<Real> LaplaceCof;
    LaplaceCof.push_back(-2 * M_PI * M_PI);
    LaplaceCof.push_back(-2 * M_PI * M_PI);
    LaplaceCof.push_back(-4 * M_PI * M_PI);

    for (int i = 0; i < numLevel; i++)
    {
        IntVect hi(IntVect::Unit * (pow(2, i + 7) - 1));
        Real dx(pow(2, -i - 4));
        Box B(lo, hi);
        FArrayBox FB(B, 2 * NumComponents);
        start = clock();
        TestExactSol2D TES2D;
        TES2D.setFAB(FB, origin, dx, 1, 0);
        TES2D.setLaplacian(FB, origin, dx, 0, NumComponents);
        end = clock();
        cout << "For sizes " << pow(2, i + 7) << "*" << pow(2, i + 7) << ", CPU time is "
             << (double)(end - start) / CLOCKS_PER_SEC << " s. " << endl;
        Box subB(lo + 2 * IntVect::Unit, hi - 2 * IntVect::Unit);
        BoxIterator subBit(subB);
        Real tempU(0), tempV(0), tempP(0);
        for (subBit.begin(); subBit.ok(); ++subBit)
        {
            IntVect iv = subBit();
            /// let numerical-laplace minus exact-laplace
            tempU += fabs(FB(iv, NumComponents) - LaplaceCof[0] * FB(iv, 0));
            tempV += fabs(FB(iv, 1 + NumComponents) - LaplaceCof[1] * FB(iv, 1));
            tempP += fabs(FB(iv, 2 + NumComponents) - LaplaceCof[2] * FB(iv, 2));
        }
        errU.push_back(dx * dx * tempU);
        errV.push_back(dx * dx * tempV);
        errP.push_back(dx * dx * tempP);
    }
    cout << " " << endl;
    cout << "3. Grid Refinement Test for 4th Discrete Laplacian:" << endl;
    cout <<"For velocity (u,v) and pressure p, we have"<<endl;
    ///Test for u.
    cout << "u:" << endl;
    for (int i = 0; i < numLevel; i++)
    {
        cout << "sizes " << pow(2, i + 7) << "*" << pow(2, i + 7) << " , L1 norm err : " << errU[i] << endl;
    }
    cout << "convergence rate is: " << endl;
    for (int i = 0; i < numLevel - 1; i++)
        cout << log(errU[i] / errU[i + 1]) / log(2) << endl;
    ///Test for v.
    cout << "v:" << endl;
    for (int i = 0; i < numLevel; i++)
    {
        cout << "sizes " << pow(2, i + 7) << "*" << pow(2, i + 7) << " , L1 norm err : " << errV[i] << endl;
    }
    cout << "convergence rate is: " << endl;
    for (int i = 0; i < numLevel - 1; i++)
        cout << log(errV[i] / errV[i + 1]) / log(2) << endl;
    ///Test for p.
    cout << "p:" << endl;
    for (int i = 0; i < numLevel; i++)
    {
        cout << "sizes " << pow(2, i + 7) << "*" << pow(2, i + 7) << " , L1 norm err : " << errP[i] << endl;
    }
    cout << "convergence rate is: " << endl;
    for (int i = 0; i < numLevel - 1; i++)
        cout << log(errP[i] / errP[i + 1]) / log(2) << endl;
}
